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Abstract 


A variable fidelity, multiscale, physics based finite element procedure for 
predicting progressive damage and failure of laminated continuous fiber rein- 
forced composites is introduced. At every integration point in a finite element 
model, progressive damage is accounted for at the lamina-level using thermody- 
namically based Schapery Theory. Separate failure criteria are applied at either 
the global-scale or the micro-scale in two different FEM models. A microme- 
chanics model, the Generalized Method of Cells, is used to evaluate failure cri- 
teria at the micro-level. The stress-strain behavior and observed failure mecha- 
nisms are compared with experimental results for both models. 


Introduction 

The ability to optimize the design of composite structures is limited by the prediction 
capabilities of numerous progressive damage and failure analysis methods. In order 
to utilize the full potential of these methods, a distinction between damage and fail- 
ure must be established. Failure indicates a global catastrophe, such as macroscopic 
cracking; therefore, it is an event that leads to large changes in the material properties 
at the failed material point. Damage, though, leads to a gradual reduction in (not 
complete absence of) load carrying capability. In most instances, damage reaches a 
critical state and becomes unstable resulting in failure; however, these mechanisms 
need to be treated separately. By considering both progressive damage and failure, 
the response of a carbon fiber laminated composite structure can be more accurately 
characterized. 

Many methods utilize failure criteria, along with linear elasticity, to predict the 
load carrying capabilities of composite structures. Typically, progressive failure is 
introduced into finite element simulations by computing load displacement behavior 
using an elasticity solution until a failure criterion is met locally at an integration 
point. Upon satisfying this criterion, specific moduli are reduced to negligible values 
at that integration point. The moduli that are abated are chosen based on the failure 
mechanism being modeled. 

This approach does not account for all the damage mechanisms observed when 
a composite structure is loaded, particularly matrix microdamage. Matrix micro- 
damage, or microcracking, in epoxy matrix composites is the growth of distributed, 
microscopic voids and fissures within the matrix, that are consequences of the man- 
ufacturing process and the deformation experienced by the polymer. Progressive mi- 
crodamage is the primary cause of nonlinearity in fiber reinforced laminates (FRLs) 
up to the onset of failure [1], As microdamage accumulates due to mechanical load- 
ing, the elastic moduli of the composite are permanently, and progressively, degraded. 
This damage manifests in the matrix of a composite structure; as a result, only E 22 
and G \2 of a lamina are significantly affected. 

Damage progression in a composite leads to a fundamentally different stress- 
strain response from the material. If the structure contains notches or holes, as is 
the case for most composite panels in service, localized reduction in moduli leads 
to redistribution of the stress and strain fields. These fields cannot be captured ac- 
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curately if a failure criterion is used in conjunction with linear elasticity because the 
properties in all locations that have not reached failure are assumed uniform. 

In a composite structure, the coupled failure of the individual constituents leads 
to the globally observed failure mechanisms. The failure of the constituents is com- 
monly modeled by reducing the lamina properties that are affected most by the prop- 
erties of the failed material. This homogenization does not capture the interaction 
between the constituents, and can lead to a misrepresentation of the failure mech- 
anisms that accrue when the structure is loaded. Micromechanics can be used to 
resolve the composite into its components. Failure can be evaluated in each of the 
constituents, thus capturing the interactions due to failure of the individual materials 
in the composite. 

Two new finite element procedures for predicting progressive damage and failure 
of FRLs were developed. The objective of this paper is to show that, for the lami- 
nates studied, the use of macroscopic failure criteria is redundant in view of assigning 
critical states to the accumulated progressive damage that occur within a lamina. 

At every integration point in each finite element model, progressive damage is 
accounted for at the lamina-level using a thermodynamically based theory developed 
by Schapery [1], Failure criteria are applied at either the global-scale or the micro- 
scale. Lamina level failure is evaluated via the 2-D Hashin-Rotem failure criterion 
[2], A micromechanics model, the Generalized Method of Cells (GMC) developed 
by Paley and Aboudi |3|, is used to evaluate the 3-D Tsai-Hill failure criterion [4] 
in the matrix phase at the micro-level, and a maximum stress criterion is utilized for 
fiber failure. Results from the two different finite element models arc examined and 
compared with experimental data from Bogert et al.[5] in Section . 


Multiscale Modeling of Damage and Failure 

Lamina Level Modeling of Progressive Damage Using Schapery 
Theory 

Thermodynamically based work potential model 

Progressive damage in the epoxy matrix composite is modeled using Schapery The- 
ory (ST) [1], This thermodynamics based, work potential theory is capable of captur- 
ing the microdamage mechanisms responsible for the material nonlinearity by divid- 
ing the total applied work, Wt into a recoverable pail (elastic), W, and a dissipated 
portion (work of structural change), IT 5 . 


W T = W + W S (1) 

As the material is loaded, a portion of the applied work facilitates structural 
changes in the material. These structural changes, such as microcracking, affect 
the elastic properties of the material. A portion of the total applied work is recov- 
ered when the structure is unloaded. The magnitude of work recovered is contingent 
upon the current degraded elastic properties. Upon subsequent reloading, the material 
will behave linearly, exhibiting the elastic properties observed during unloading, until 
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the material reaches the previous maximum strain state. After this state is achieved, 
structural changes resume, further degrading the elastic moduli of the material. This 
process is shown in Figure 1. The shaded area represents Wg, and the area under 



Figure 1 . Irrecoverable portion of total work, Ws, is represented by the shaded area. 


the linear unloading curve is W. It is assumed that the material behaves as a secant 
material, a reasonable assumption for FRLs [6], 

Both W and Ws are functions of a set of internal state variables (ISYs), S m , (to = 
1,2, M). These ISVs account for any inelastic structural changes in the material. 
Differentiating Ws with respect to any ISV, S m yields the thermodynamic force, f m , 
available for producing structural changes associated with the m th ISV. 


_ dW s 
dS m 


(2) 


It is shown in Ref. [1] that the total work is at a minimum with respect to each ISV. 


dW T 

dS, „ 


(3) 


Additionally, Rice [7] has shown that according to the second law of thermodynam- 
ics: 

frnSra > 0 (4) 

Equations (2), (3), and (4) form the foundation of a thermodynamically based work 
potential model for nonlinear structural changes in a material. 


Application of ST to fiber reinforced plastic composites 

Damage accumulates in a composite through numerous mechanisms affecting the 
constituent materials in the composite. Micro-level damage, which includes matrix 
microcracking, fiber kinking and debonding, is a class of damage mechanisms sep- 
arate from matrix failure due to transverse (macroscopic) cracking. Matrix damage 
accumulates gradually at the micro-level until its effects are superceded by the rapid 
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progression of macroscopic cracks. Like matrix macrocracking, fiber damage also 
does not typically occur progressively, but rather, abruptly. Moreover, once the fibers 
in a composite lamina begin to break, it has nearly lost its entire load carrying capa- 
bilities. Across the breakline, however, the adjacent layers carry the loads through 
load re-distribution. Unfortunately, microdamage mechanisms are often overlooked 
in analyses; instead, the more catastrophic matrix macrocracking and fiber breakage 
are the focus. ST is capable of modeling the effects of progressive microdamage in 
the matrix phase of FRLs. 

The inelastic work of structural change, Wg, can be a function of any number of 
state variables. To apply this work potential model to progressive failure analysis, 
it is assumed that the structural changes which result from microdamage, depend on 
only one ISV, S. This ISV is assumed responsible for all material nonlinearities up to 
macroscopic matrix cracking (global cracking), delamination and fiber breakage, and 
account for all damage present in the matrix of each lamina in a composite structure. 

It can be assumed that Wg is an additive function of the ISVs, Wg = J2T Ilu ( .5); ) . 
Furthermore IT) are in one-to-one correspondence with their arguments; so, if) can 
be chosen such that Wi = S). Since, in this case, Wg is a function of only one ISV, 
Wg can be chosen such that ISV, Wg = S. Therefore, the ISV governing the amount 
of work used to advance microdamage and that actual work are equivalent. Equation 
(1) can now be recast. 

W T = W + S (5) 

Differentiating (5) with respect to S, and utilizing Equation (3) yields: 


dW 

~dS 


( 6 ) 


Additionally, combining Equations (2) and (4) with Wg = S results in 


S > 0 


(7) 


which is a statement on the inadmissibility of damage “healing”. Equation (7) dictates 
that the amount of work used to progress microdamage can never decrease; therefore, 
that energy has been dissipated into creating structural change and cannot be recov- 
ered. The combination of Equations (6) and (7) represent the evolution equations for 
microdamage in the matrix of the composite. 


Formulation of constitutive law 

Plane stress lamina stress-strain relationships can be written in principal material co- 
ordinates as 

Oil = Q llCll + Q 12^22 

022 = Q 12^11 + <^22^22 (8) 

T 12 = <366712 
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where 712 is the engineering shear strain and 


Qu — 

Q22 = 

Q 12 = 
^66 = 

U 2 l = 


E 11 

1 -V 12 V21 

E 2 2 

1 - ^ 12^21 
^12^22 
Gi 2 
V12E22 
E u 


( 9 ) 


where Efi 1 is the axial elastic modulus, E22 is the transverse elastic modulus, 77 2 is 
the Poisson’s Ratio, z/ 2 i is the transverse Poisson’s Ratio and G , 9 is the elastic shear 
modulus. After assuming that the quantity ^ 12 z/ 2 i << 1 , Equations ( 9 ) become, 


Qn = E 11 

Q22 = E22 

Q\2 — V12Q22 

Qg 6 ~ Gl 2 


( 10 ) 


Determining the damage state 

It is necessary to define in what manner the moduli degrade as functions of the ISV. 
Since the damage mechanisms considered in this progressive damage approach are 
exclusive to the matrix of the composite, it is safe to assume that the moduli affected 
by this damage are limited to E 2 2 and G i 2 . These moduli can be written as functions 
of S. 

E 22 = E22oe s (S ) ( 1 1 ) 

G 1 2 = Gi2og s (S ) ( 12 ) 

where /f 220 and G i 20 are the undamaged transverse and shear elastic moduli, e s (S) 
and g s (S ) are factors relating the transverse and shear moduli to the microdamage, 
S. Sicking [6] provided a procedure for determining e s and g s experimentally. The 
experimental curves can then be fit with polynomials (such that moduli at S' = 0 
arc E22 = E220 and G 12 = G120. corresponding to an undamaged state) and used in 
Equations ( 1 1 ) and ( 12 ). 

The elastic strain energy density, W, can be written using the constitutive rela- 
tionships. 

II = -(E n €^ + £’ 22 e 22 + Gi 2 7 ^ 2 ) + Qi 2 ene 22 (13) 

Employing Equation (6) with ( 1 1 ), ( 12 ), and ( 13 ), and assuming the quantity Q12 = 
^i2f?22 is constant and independent of S, yields a damage evolution ordinary differ- 
ential equation which can be solved for S. 

e 22 ®E 2 2 7 12 dGl 2 = , 

2 dS 2 8 S { ’ 

The above equation indicates that the work of structural change depends only on 
the strain state, the initial virgin composite moduli (f? 22 0 and G i 20 ), and the damage 
functions (e s and g s ). 
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Experimentally it has been determined that S behaves as e 3 , thus it is convenient 
to introduce a reduced damage variable, S r . 


S r = S 1/3 


(15) 


Using the reduced ISV the evolution equation, Equation (14), becomes 


e| 2 dE 2 2 7i2 dG 12 

2 dS r 2 8S r 


(16) 


Once S is determined from Equation (16), the transverse and shear moduli can be 
degraded accordingly using Equations (11) and (12). 

In previous studies [8], where the loading was compression dominated, the in- 
stantaneous fiber rotation is captured in conjunction with ST. This has the advan- 
tage of predicting fiber kinking failure, avoiding the use of an explicit fiber direction 
compressive strength criterion (see Equation (45) later). In the present work, where 
loading is considered tension dominated, the effect of including the fiber rotation was 
initially deemed unnecessary. However, the discoveries reported in Section indicate 
that including this feature would also be worthwhile for this problem. 


Micromechanical Modeling Using the Generalized Method of Cells 

A micromechanical analysis technique, coined the Method of Cells, was developed 
by Aboudi [9]; subsequently, Paley and Aboudi [3] expanded the Method of Cells 
into the Generalized Method of Cells (GMC), and later Aboudi et al. [10] fur- 
ther increased the robustness of this method with the High Fidelity Method of Cells 
(HFGMC). These methods provide semi-closed form solutions for determining global 
anisotropic composite properties in terms of the constituent materials, as well as, 
stresses and strains in each of the constituent subcells. The sophisticated methods 
(GMC and HFGMC) offer a high degree of accuracy at a relatively low computational 
cost. The following sections detail the formulation of GMC (employed herein). The 
reader is referred to Ref. [10] for details on HFGMC. 

Kinematics and constitutive relationships 

It is assumed that a unidirectional fiber composite can be represented as a collection 
of repeating volume elements (RVE). Paley and Aboudi[3] chose to model this RVE 
as an element consisting of Ng x N y (/), 7 = 1, 2, Ng r ,) subcells as shown in Figure 
2. Each of these subcells is occupied by one of the constituents in the composite. The 
number of subcells and the materials occupying each subcell is completely general. 
For a two-phase fibrous composite any desired micro-structure can be represented by 
occupying each subcell with either a matrix or fiber constituent. 

The aq-axis shown in Figure 2 is the fiber direction, and the cross-sectional area 
of each subcell is given by hgt 7 . A local coordinate system can be 

introduced with its origin located at the center of each subcell, as shown in Figure 3. 
The objective of this method is to determine the average behavior of the composite 
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Figure 2. Representative volume element used in GMC [3], 



x 3 


Figure 3. Local coordinates used in GMC subcells [3], 
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material; thus, it is sufficient to model the displacements in each subcell using a linear 
theory (HFGMC employs a higher order displacement approximation). 


uf 7) = wf Y) + xfUT y) + PP4 0l) 


( 17 ) 


where i = 1,2,3, w. 
ables 


( 07 ) 


is the displacement at the center of subcell f3^. Microvari- 


' ) characterize the first-order dependence of the displacement field 
on the local coordinates xp and xp. 

The components of the strain tensor follow from Equation (17) as, 


,(07) 

t ij 


( diuf 7) + djuf 1 ^ 


( 18 ) 


where d\ = p - , d 2 — Swr ■ — — rvr- Substituting Equation (17) into Equation (18) 

ax l ctej cteg 

results in the six components of the average strain tensor for each subcell in terms of 
the microvariables. 



diw[ 0 ^ 




efP = ^ 

2e-g 7) = ^ + 4" 7) 
2efp 1 = p 0l) + d lW f ' y) 


( 19 ) 


2e[P = p 0l) + d lW f l] 

The constitutive law also needs to be defined for each subcell, and can use any 
stress-strain relationship desired. In this work, it will be assumed that only elastic 
subcell strains, ppx are present. However, the constitutive law can be amended to 
incorporate additional strains such as thermal and inelastic strains. Hooke’s law for 
relating the subcell stresses, ap : } , to the elastic subcell strains can be written as 


=.( 07 ) _ fiiPl) =( 07 ) 
°ij ~ L ' ijkl t kl 


where C, 


( 07 ) 

ijkl 


are the components of the elastic stiffness tensor. 


( 20 ) 


Displacement continuity conditions 

It is required that subcell displacements are continuous at the interfaces between ad- 
jacent subcells, as well as at the boundaries between neighboring repeating cells. 
Enforcing these conditions will yield 2 (Ng + iV 7 ) + NgN^ + 1 equations. Moreover, 
all microvariables are eliminated from these equations. For detailed derivations of 
displacement and traction continuity conditions see Ref. [3], 
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The first necessary equation is obtained by defining the average strains in the 
composite, , in terms of the average subcell strains. 


n b n 7 


e ij ~ i p X/ X] 

p=\ 7=1 


Mi) 


( 21 ) 


where h, £ and hp, L, represent the repeating cell, and subcell geometry, respectively 
(see Figures 2 and 3). 

Next, displacement continuity is satisfied in an average sense over all subcell and 
repeating cell interfaces. Since these continuity conditions are satisfied on average, 
the shape of the fiber does not appeal - in the final result. Thus, no stress concentra- 
tions are developed at the corners, and the end result is that the subcell strains, and 
stresses, are determined as a function of only the fiber volume fraction and constituent 
properties. 

After enforcing average displacement continuity on the interfaces, the following 
2 (Np + iV 7 ) equations are produced: 


Nf 3 

J2 hpe = he-22, 

p=i 

iv 7 

Vs3 7) = ^33, 

7=1 

Np 

^/3^i2 7 ^ = hei2, 

P=1 




y, ^7 e 13 7) — ^ e 13, 

7=1 


7 = 1,...,-ZV 7 

P — 1 ) — -5 Np 

7 = l,...,iV 7 
Np 


( 22 ) 

(23) 

(24) 

(25) 


The final A 7 A 7 equations come from enforcing uniform strain in the xj -direction 
across all subcells. 


where 

and 


es = 


£11 


(26) 

md rewritten in matrix from. 


J 6 


(27) 

£23, 2 ei 3 , 2612 } 

(28) 

. . . , e {N e N ^] 


(29) 


where are vectors containing the average subcell strains in the same order as 
Equation (28). 
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Traction continuity conditions 

Traction continuity must also be enforced in order to arrive at the correct number of 
equations needed to solve for the 6 NgN^ subcell strain unknowns. However, some of 
the traction continuity conditions are redundant. After eliminating repeating traction 
conditions, the following §NgN 1 — 2 (Ng + AT) — 1 independent equations remain: 


+ 77 ) 

°22 

- 7 (/37) 

— °22 > 

P = 

1,- 

..,Ng- 

• 1, 7 = 

1 ,...,a 7 

(30) 


= 

P = 

1,- 

■■,Ng, 

7 = 1,.. 

• • A- - 1 

(31) 

+ 7) 

= 47, 

P = 

1,- 

■ ■ j Ng — 

1, 7 = 

1,...,1V 7 

(32) 

+/A) 

°32 

- 7 (/37) 
“ a 32 , 

P = 

1,- 

■■ ,Ng . , 

7 = 1,.. 

,,1V 7 -1 

(33) 

+/A) 

a 21 

~ a 21 ) 

P = 

1,- 

..,Ng- 

1, 7 = 

1 ,...,a 7 

(34) 

+31 

- 7 (/37) 
“ +31 ? 

P = 

1,- 

■■,Ng, 

7 = 1,.. 

,,1V 7 -1 

(35) 


where (3 and 7 are given by 


P = 


P + 1, P < Ng 
1, P = Ng 


(36) 


7 = 


(37) 


7 + 1, 7 < iV 7 
{ 1, 7 = ^7 

These traction conditions can be recast in terms of the average subcell strains 
using the constitutive relationship, in this case Equation (20). The equations can then 
be rewritten in matrix form. 


A Meg — 0 


(38) 


where eg is given in Equation (29) 


Determining subcell strains 

Once Aq, A m , and J have been determined, the subcell strains can be computed by 


solving 





Ae s 

= Ke 

(39) 

where 





A = 

Am 
. A g . 

(40) 

and 





K = 

0 

J 

(41) 


After the subcell strains are obtained, it is trivial to produce the subcell stresses using 
the constitutive law. 
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Failure Criteria 

Modeling lamina-level failure with Hashin-Rotem failure criteria 

In this work, failure is treated separately from damage. Damage is considered the pro- 
gressive deterioration of the matrix. A method for modeling damage was presented in 
Section . Failure, however, is a localized, catastrophic event, after which the material 
can no longer support any load in a particular direction at the point of failure. 

There are numerous theories in the literature that offer methods for evaluating 
the onset of catastrophic failure. One criterion used in this investigation is the 2-D 
Hashin-Rotem (H-R) failure theory [2], which is also used in Ref. [5], This lamina 
theory is the simplest, empirical theory that still retains the transversely isotropic 
nature of unidirectional fiber composites. 

Four separate criteria encompass the H-R failure theory. Transverse (perpendic- 
ular to the fiber direction) and shear stresses are responsible for matrix failure in a 
unidirectional fiber composite. Matrix failure at any material point is dictated by 

(^) 2 +(^V)‘ = <£, ff22 > o (42) 

when the transverse stresses, cr- n , are greater than zero (tension), and 

(f ) 2 + (?) 2 = *- ^ <0 (43) 

when (T 22 is less than zero (compression), where t \2 is the shear stress in the composite 
material coordinates, Y t is the matrix transverse strength in tension, Y c is the matrix 
transverse strength in compression, and T is the matrix shear strength. Once d m is 
greater than or equal to one, that matrix at that location has failed, and that point can 
no longer sustain transverse or shear loads. Although, it is assumed that such a point 
will still be able to carry loads in the fiber direction until it reaches the limit of fiber 
failure. 

A similar set of criteria is used to govern fiber failure. When the stress in the fiber 
direction, a n , is greater than zero, the fiber failure criterion used is 

(j^y = d f’ °n>0 (44) 

Whereas, when the fibers in the composite are subjected to a compressive load (<Tn < 
0) 

Gr) =d *’ an<0 (45) 

is used. The fiber direction tensile strength is X t , and X c is the compressive strength. 
Fiber failure occurs when d f > !? after which, that point can no longer carry axial 
loads. 

The H-R failure criterion provides a simple 2-D model for predicting localized 
failure in a finite element model. However, as is mentioned in Ref. [5], the predictions 
of the failure criteria are dependent on the density of the finite element mesh near a 
notch tip. This indicates that new stress measures that incorporate the length scales 
associated with the discretization process needs to be developed in order to obtain a 
robust prediction model [11]. 
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Modeling micro-level constituent failure with Tsai-Hill 


The multiscale computational method (MCM) employed in one of the FEM sim- 
ulations evaluated uses GMC to resolve the applied global stresses to the micro- 
constituent level. Failure of the matrix phase is established using the 3-D Tsai-Hill 
failure criterion[4]. Assuming failure in the matrix subcells is isotropic, the criterion 
is given as 


(ff<f ” 7 ")) 2 + (ct ^” 7 ” 1 ) 2 + (ctS ” 7 " 1 ) 2 


Ylt 

XPm'Ym) — (Pm'Ym) 

"°11 °22 


<7 


(/3m7fO - (/?m7m) - (/^m7m) 


11 


a- 


33 


a. 


22 


a- 


33 


V 2 

± mt 

( -z.(0m7m) \2 I /Adm7m)\2 i \ 2 

l°12 ) t°13 ) (°23 ) 


j-*2 


— <722 > 0 


(46) 


where / 3 m and 7 m are the matrix subcell indices and the average applied transverse 
stress, < 7 2 2 , is tensile. Similarly, for applied compressive transverse stresses: 


(-gm7m))2 + ^iPrnlm)^ + 

Y 2 

me 

— a 11 a 22 — Oil ^33 — a 22 °33 


Y 2 

me 

(-gm7m))2 + (^7m)j2 + (^g-7m))2 

Y~2 

L m 


— dm, °22 < 0 


(47) 


where Y mt and Y mc are the matrix transverse tensile and compressive strengths, re- 
spectively, and T m is the matrix shear strength. 

Maximum stress criteria , analogous to Equations (44) and (45), are used to dictate 
failure in the fiber subcells. 


JfinfY 

°ii 


X 


ft 


= d% 


(XU > 0 


(48) 


/=.(d/7/)\ 2 

(^r) =4. *.<0 (49) 

where ffj and 7 f are the fiber subcell indices, Xj t and Xf c are the tensile and com- 
pressive fiber strengths, and d\ \ is the applied axial stress. 

Failure arises in a matrix subcell when d m > 0 and in a fiber subcell when d f > 
0 for that subcell. Once a subcell has failed, all the properties of that subcell are 
degraded appropriately. 


Finite Element Model Description 

Two different finite element models were used to predict the behavior of a notched 
carbon fiber reinforced epoxy panel first introduced in Ref. [5], Two laminate stack- 
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ID 

Stacking Sequence 

Thickness (in.) 

Laminate Sequence- 1 

[0] 12 

0.078 

Laminate Sequence-2 

[45/0/-45/0/90]s 

0.065 


TABLE I. Laminate stacking sequences investigated. 


Property 

Value 

E u (Msi) 

23.2 

E 2 2 (Msi) 

1.3 

G \2 (Msi) 

0.9 


0.28 


TABLE II. Initial elastic properties of T800/3900-2 lamina. 

ing sequences, shown in Table I, were modeled, and the elastic properties correspond- 
ing to T800/3900-2 (Table II) were used as the initial properties for each layer. 

The mesh used in these FEM models, shown in Figure 4, consists of 6642 nodes 
and 6480 Abaqus S4R elements. This mesh was used because the high density of ele- 
ments near the notch tips is needed to produce accurate stress fields at those locations 
[ 12 ]. 

No constraints were placed on the vertical edges of the model, and the bottom 
edge was restricted from moving in the x, y, z, and all rotational degrees of freedom. 
The top edge is fixed in the x, and z displacements and all rotations. A vertical 
displacement is applied in the y direction to simulate tensile loading. 

Static analysis is performed in Abaqus/Standard, and an edge displacement of 
0.025 inches and 0.065 inches is applied to laminate stacking sequence 1 and 2, re- 
spectively. The maximum allowable displacement in each time step was set to 0.1%, 
and the minimum allowed displacement was 0.1E-7%. 

The coupled microdamage-failure models were implemented using the Abaqus 
user subroutine, UMAT [13]. Both models use strain to calculate the reduced damage 
state, S r , with Equation (16). The damage functions, e s and g s given in Equations 




(a) FEM Mesh used in simulations. 


(b) Enlarged view of mesh near notch tip 


Figure 4 
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E22 coefficients 

Values 

G 1 9 coefficients 

Values 


1.0000 

9 s 0 

1.0000 

&sl 

-0.035 1 

9 si 

-0.0377 

%2 

-0.0096 

9s2 

-0.0237 

e s3 

-0.0016 

9 s 3 

0.0053 

%4 

0.0003 

9sa 

-0.0004 


TABLE III. Microdamage polynomial coefficients for E 2 2 and G 12 . 


Property 

Value (Msi) 

Y t 

0.00872 

Y c 

0.0243 

T 

0.0048 

x t 

0.412 

X c 

0.225 


TABLE IV. Hashin-Rotem failure strengths. 


(11) and (12) were modeled as fourth order polynomials. 

e s — e so + e si S r + e S 2S~ + e S 3 + eg^Sr (50) 

9 s — 9s0 + 9 s 1 Sr + 9 32^1 + g s ?,Sr + 9s4,S^ (51 ) 

The polynomial coefficients in Equations (50) and (51) are given in Table III. These 
values were obtained by scaling the values reported in Ref. [6] by the ratio of the 
respective virgin elastic moduli. 

It is possible that the ratio of G 12 /E 22 , as calculated by Equations (11) and (12), 
reaches a value that produces an ill-conditioned Jacobian in the FEM solution. There- 
fore, G 12 is restricted from falling below 40% of G V2 0 , which corresponds to ,5' r =3.74 
psis . After any point reaches this state, S r is still calculated but it is not used to update 
G 12 nor E 2 2 - Exploration into methods to avoid this instability and still allow dam- 
age to progress are currently underway. However, the possibility that this numerical 
instability is signaling a critical damaged state is also being investigated. 

Method 1 - Schapery Theory with Lamina Level Hashin-Rotem 
(ST/H-R) 

The first algorithm employs the H-R failure criteria at the lamina level. At each 
material point, the failure criteria are evaluated using the properties in Table IV. All 
failure strengths were taken directly from Ref. [5], except the shear strength, T. The 
shear strength was adjusted to account for the shear modulus degradation resulting 
from the progressive damage in the laminate. The value that yielded the best results 
for laminate sequence- 1 was used for both laminates. 

If Equation (42) or (43) is satisfied, E 2 2 , G 12 , and v X2 are reduced to 35% of 
their current values. Each subsequent time step, the calculations involving ST are 
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Figure 5. 2 x 2 RYE used in micromechanics simulations. 


Fiber Properties 

Values 

Matrix Properties 

Values 

E^ (Msi) 

42.49 

E™ (Msi) 

0.3415 

E{ 2 (Msi) 

13.2 

E% (Msi) 

0.3415 

u 12 

0.2316 


0.35 

—J 

V 21 

0.45 


0.35 

G{ 2 (Msi) 

8.0 

G™ 2 (Msi) 

0.3267 


TABLE V. Elastic properties of fiber and matrix constituents used in GMC. 

circumvented and E 22 , G 12 and u 12 are further reduced by 35% until E 2 2 and G 12 fall 
below 1000 psi, and is 12 is less than 0.001. 

After matrix failure has occurred, Equations (44) and (45) are still evaluated to 
check for fiber failure. Once either of those criteria are satisfied, E n is also reduced 
to 35% then reduced by the same percentage in each following time step until it falls 
below 10,000 psi. 

Method 2 - Schapery Theory with Micro-level Tsai-Hill (ST/GMC/T- 
H) 

The second method for modeling damage coupled with failure follows the same pro- 
cedure for determining the damage in the structure as the previous method; however, 
the failure criteria is no longer evaluated at the lamina level. Instead, the MAC/GMC 
suite of micro-mechanics codes is used [14, 15], At each material point, MAC/GMC 
is called and a 2 x 2 RVE, as shown in Figure 5, is used to model that point. The 
RVE consists of three matrix subcells and one fiber subcell. The fiber and matrix 
constituents have the initial properties given in Table V. To provide global proper- 
ties consistent with those given in Table II anisotropic properties were used for the 
constituent cells. 

At every time step, the micromechanics model must produce composite moduli 
that are consistent with the lamina level moduli calculated using Equations (IT) and 
(12). Therefore, it is necessary that the moduli of the matrix constituents degrade in a 
manner that produces consistent E 22 , and G ± 2 values. Two fourth order polynomials 
are used to calculate the matrix Young’s Modulus and shear modulus, E m and G m , as 
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E m coefficients 

Values 

G m coefficients 

Values 

CmO 

1.0000 

9m.O 

1.0000 

Cm. 1 

-0.0378 

Qml 

-0.0918 

Cm.2 

-0.0128 

9 m2 

-0.0560 

Cm3 

-0.0005 

9 m3 

0.0137 

^m4 

0.0002 

9mA 

- 0.001 1 


TABLE VI. Microdamage polynomial coefficients for E 22 and CI- 


Property 

Value (Msi) 

Ymt 

0.007 

Y 

1 me 

0.018 

T rn 

0.006 

Xft 

1.100 

X .fc 

0.414 


TABLE VII. Constituent failure strengths. 


a function of S r . 


E m — E m o(e m Q + e m i S r + e m 2 + C m ^S^ + e m /±S ' r ) 


(52) 


iri E m () i fjrnt) V fjrri 1 E r V fjrn'2'^r % 9m3^ r "E 9mlE r ) 


(53) 


where E m0 and G m0 are the undamaged matrix stiffnesses, and the polynomial coef- 
ficients are given in Table VI. The matrix Poisson’s Ratio, v m , remains unchanged. 


Using consistent properties for the matrix modulus, the subcell stresses are cal- 
culated. Using these stresses, micro-level failure criteria are evaluated. Since the 
micromechanics model is 3-D, the 2-D H-R failure criteria are no longer applicable. 
Instead, 3-D Tsai-Hill (T-H) failure criteria are used for the matrix subcells 

The constants used in Equations (46)-(49) are given in Table VII. The matrix 
shear strength, T m , was obtained by applying a global shear stress, 072 , to the RVE. 
The resulting largest matrix subcell shear stress, a 12 is used as the matrix shear 
strength. Similarly, Y mt and Y mc are determined by applying global transverse stresses 
a 22 = Y f and Y c . Equations (46) and (47) are then solved with d rn = 1, 6 = 1, and 
7 = 2 for the matrix axial strengths in tension and compression. Finally, Xf t and Xj r 
are obtained by subjecting the RVE to global axial stresses <7n = X t and X c . The 
ensuing axial stresses in the fiber subcell, d^ 1 ' ) are used for maximum allowable fiber 
stresses in tension and compression. 

If any of Equations (46), (47), (48), or (49) are met in any matrix subcell, all the 
subcell properties are reduced by 99%. As in method 1, once failure has occurred in 
any subcell, progressive damage is deactivated at that integration point. 
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Laminate 

Load Type 

Experiment 

ST/H-R 

ST/GMC/T-H 

Laminate 1 

Splitting 

8635 lbs. 

8600 lbs. 

8792 lbs. 

Laminate 2 

Ultimate 

12632 lbs. 

12740 lbs. 

12322 lbs. 


TABLE VIII. Critical Loads. 

Results and Discussion 

Load versus edge displacement data is reported for all three laminates using both 
models (ST/H-R and ST/GMC/T-H) and compared to experimental results. Addi- 
tionally, applied local strains are compared to experimental results obtained from 
four strain gages placed at the locations shown in Figure 6(a). The elements used to 
represent these strain gages are shown in Figures 6(b) and 6(c). Data in these ele- 
ments are averaged to obtain the strain data at that strain gage location. Table VIII 
contains landmark experimental and computational loads in each of the laminates, 
such as splitting and ultimate loads. 



(a) Strain gage locations (b) Elements used to represent strain gages in simulation, 

in experiment [5]. 



(c) Elements used to represent 
Sg-4 in simulation. 


Figure 6 


NASA/TM- 


-2008-215448 


18 



Figure 7. Load versus displacement of a 4” section, laminate sequence 1 . 


Laminate Sequence- 1 

Applied load versus the edge displacement (of a 4 inch section) results from both 
techniques are compared to experimental data in Figure 7. It can be seen that the 
multiscale method, ST/GMC/T-H, yields data that is extremely close to experiment. 
The load/displacement response obtained using ST/H-R are still desirable. Although 
the curve exhibits softening not present in the experiment or in the other model, this 
softening is well after the experimental splitting load, which was reported to be 8635 
lbs. in Ref. [5], and may be a product of the semi-gradual manner in which the 
properties of the failed elements are reduced. 

Far field load versus local strain data is presented for all four gages in Figure 

8. The splitting load predicted by both simulations can be extracted from the sg-1 
data (Figure 8(a)). ST/GMC/T-H predicts a splitting load of 8600 lbs., and ST/H- 
R anticipates a splitting load equal to 8792 lbs. These values are reported in Table 
VIII. Localized hardening can be observed in the results. This hardening is actually 
localized strain relaxation at sg-1 and is a consequence of failure near the notch. The 
data at sg-4 using both approaches also displays some variation from the experimental 
data. This is to be expected, however, because the strain gradients near the notch tips 
are extremely high. 

The failure patterns produced by ST/GMC/T-H and ST/H-R are shown in Figure 

9. Matrix failure is represented in the ST/H-R results with green elements, and fiber 
failure is indicated with red elements. The ST/GMC/T-H results represent 1 matrix 
subcell failure with turquoise elements, 2 failed matrix subcells with green elements, 
3 failed matrix subcells with yellow elements, 1 failed fiber subcell and 3 failed matrix 
subcells (4 total) with red elements. 

The matrix failure patterns generated from both simulations are in agreement with 
each other and represent matrix shear splitting exhibited in experiment (Figure 10). 
The two methods produce different fiber failure patterns. At the onset of fiber fail- 
ure, stiffness contributions due to the matrix near the notch are severely diminished. 
However, it is possible that some of the matrix in some elements near the notch is 
still intact. Additionally, initial fiber failure can further contribute to the variations 
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P (Ibf.) P (Ibf.) 




(a)sg-l (b) sg-2 




(c) sg-3 (d) sg-4 


Figure 8. Load vs. strain for Laminate 1. 


NASA/TM— 2008-2 1 5448 


20 


(a) ST/GMC/T-H failure pattern, P= 14048 
lbs. 


(b) ST/H-R failure pattern, P= 14082 lbs. 



(c) ST/GMC/T-H damage pattern. P = 8600 (d) ST/H-R damage pattern. P=8792 lbs. 

lbs. 


Figure 9. Failure and damage paths in 0° layer of Laminate 1 . 


Notch 



Damage Path (Splitting) 

Figure 10. C-Scan of failed laminate 1 specimen exhibiting splitting [5], 
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(a) ST/GMC/T-H 


(b) ST/H-R 


Figure 1 1 . Magnification of failure in laminate 1 at notch tip. 



Figure 1 2. Load versus displacement of a 4” section, laminate sequence 2. 

in stiffnesses near the notch tip. This gradient in material properties can lead the 
notch tip to enter a bifurcated state (Figure 1 1(b)), which in turn drastically affects 
the stress distribution in the surrounding region. This a facet of the simulation and 
is not occurring in the experiment because the material is actual opening at the notch 
tip and in the split regions. Additionally this behavior is highly sensitive to the ratio 
of the matrix strengths, and displays the highly coupled nature of failure mechanisms 
in composite materials. Thus, the fiber failure displayed by ST/H-R is not physical. 
The fiber failure path produced by ST/GMC/T-H is the expected mechanism (Figure 
11(a)). 

Figure 9 also shows the microdamage patterns calculated using both methods. 
These damage patterns are similar, and they progress in the same manner as the ma- 
trix failure shown in Figures 9(b) and 9(a). This indicates that matrix failure in 0° 
laminates is a culmination of microdamage. 

Laminate Sequence-2 

Figure 12 shows the bulk response for laminate sequence-2 and Figure 13 show the 
local strain gage data for both simulations and experiment. The ultimate load predic- 


NASA/TM— 2008-2 1 5448 


22 


P (Ibf.) P (Ibf.) 




(a)sg-l (b) sg-2 




(c) sg-3 


(d) sg-4 


Figure 13. Load vs. strain for Laminate 2. 
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(c) ST/GMC/T-H damage pattern. P = 4878 (d) ST/H-R damage pattern, P=5572 lbs. 

lbs. 


Figure 14. Failure and damage paths in 45° layer of Laminate 2. 

tion for both ST/H-R (12322 lbs.) and ST/GMC/T-H (12740 lbs.) are very close to 
the reported experimental ultimate load (12632 lbs.). Again, though, there is a degree 
of softening in the simulations not present in the experiment. As the fibers fail near 
the notch in the model, there is no stiffness in any direction at those points, which is 
representative of a crack propagating out from the notch tip. In the experiment, this 
crack initiates at the notch tip and propagates rapidly to the edge of the specimen. 
Since all FEM calculations are preformed using an implicit solver, it is not possible 
to capture the rapid progression of the crack. Therefore, its is the gradual nature of 
the fiber failure and progression of the crack tip that yield the exhibited softening in 
the models. This is also evident in the local strain gage results. The initial nonlin- 
earity observed in the experiment is captured; however after premature fiber failure 
initiates, the local stress-strain curves begin to deviate 

The failure paths for the +45° layer are displayed in Figure 14 and use the same 
convention as the previous laminate. The failure patterns for the top layer closely rep- 
resent the failure mechanisms observed in experiment (Figure 15). However, there is 
failure mechanism present in the ST/H-R simulation not present in the multiscale 
calculations. The laminate modeled using ST/H-R experienced fiber failure, due to 
compression at the point where the notch radius meets the slot, see Figure 1 6. Figure 
17 shows the failure patterns in the 0° and 90° plies. The presence of compressive 
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Figure 15. Photograph of failed Laminate 2 specimen [5]. 



Figure 16. Compressive fiber failure in 45° layer of laminate 2. 



(a) ST/GMC/T-H 
failure pattern in 
0° layer. 


(b) ST/H-R fail- 
ure pattern in 0° 
layer. 


(c) ST/GMC/T-H 
failure pattern in 
90° layer. 


(d) ST/H-R fail- 
ure pattern in 90° 
layer. 






(e) ST/GMC/T-H 
damage pattern in 
0° layer. 


(f) ST/H-R dam- 
age pattern in 0° 
layer. 


(g) ST/GMC/T-H 
damage pattern in 
90° layer. 


(h) ST/H-R dam- 
age pattern in 90° 
layer. 


Figure 17. Failure and damage paths in 0° and 90° layers of Laminate 2. 


NASA/TM— 2008-2 1 5448 


25 








fiber failure in this laminate demonstrates that incorporating fiber rotation[8] in the 
progressive damage calculations is warranted, even though the specimens are loaded 
in tension. An added advantage of this would be the ability to dispose of a failure 
criterion to predict compressive fiber failure by setting a limit on the amount of al- 
lowable fiber rotation or rate of fiber rotation. Both simulations exhibit the same 
matrix splitting and fiber failure in the 0° plies and an extensive amount of matrix 
damage in the 90° plies. 

Microdamage patterns, produced far before the accumulation of failed elements 
is significant, show patterns similar to the matrix failure patterns produced when the 
laminate has reached its ultimate load. The microdamage patterns for all layers are 
given in Figures 14 and 17. This indicates that matrix failure in this laminate is a 
product of microdamage and may be characterized by incorporating a critical limit 
on the ISV, S r , rather than using a macroscopic failure criterion for the matrix. 

The use of constant failure strengths may be influencing the results of laminate 
2. Matrix shear strength, T, governs the failure behavior in the 0° lamina. Since this 
layer carries more load than any other ply in the the multi-angle laminate, the ultimate 
loads of the laminate is dictated by the performance of this lamina. However, T and 
T m are calibrated using laminate sequence- 1 . At the onset of matrix failure, the failed 
elements have already damaged progressively; hence, the material at the location of 
matrix failure is different than the material in the far field. Furthermore, the level of 
progressive damage near the notch in all layers of laminate 2 is different than that 
experienced in laminate 1 ; therefore, the material properties at those locations will 
also be different. Using the same shear strength for all laminates is analogous to 
predicting failure in different materials with the same shear failure constant. This 
argument can be extended to the other failure constants influencing matrix failure 
( YtXcXmt , and Y mc ): although the influence of these parameters on the bulk response 
of the laminates is not as profound as T or T m . 

Test results for a third laminate sequence were reported in Bogert et al. [5]. How- 
ever, in this third sequence, significant delaminations were observed, and it was re- 
ported that progressive failure in the laminate occurred simultaneously with delami- 
nation failure. In the present work delamination as a failure mechanism is not mod- 
eled. A model for delamination, however, is presented in Ref. [16], 


Conclusion 

Progressive matrix microdamage in a uniaxially loaded center notched specimen has 
been modeled using Schapery Theory (ST). However, the damage modes associated 
with progressive damage are separate from those caused by catastrophic failure. Thus, 
it is necessary to develop a criterion that evolves progressive damage into catastrophic 
failure. Two novel methods for capturing both progressive damage and failure were 
presented. One uses lamina level ST and lamina level Hashin-Rotem (H-R) failure 
theory (ST/H-R). The other incorporates lamina level ST with failure evaluated at the 
micro-constituent level using the Tsai-Hill (T-H) failure criterion for the matrix and 
a maximum stress criterion for the fiber (ST/GMC/T-H). Failure strengths used in 
these criteria were calibrated against experimental results from a [0] 12 T800/3800-2 
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laminate. 

The performance of these methods were then evaluated on a multi-angle laminate 
([+45/0/-45/0/90]s) of the same material. Additionally the progressive damage pat- 
terns mirrored the matrix failure patterns determined using a matrix failure criterion. 
Thus, this criterion may be redundant, and it may be possible to eliminate the ma- 
trix failure criterion altogether by determining a critical damage state, evolved from 
the progressive ST damage model. Additionally, fiber failure in compression is con- 
trolled by fiber rotation, which is a product of the stiffness in the surrounding matrix. 
Therefore, this failure mechanism can also be predicted using progressive damage, 
as shown earlier in Ref. [8], Since, fiber failure in tension is not a progressive phe- 
nomenon, a method to account for fiber breakage in tension is still needed. There- 
fore, this study has shown that accurate modeling of the progressive degradation of 
the matrix, together with a method to account for fiber direction failure in tension and 
a critical evolved damage state are sufficient to completely capture in-plane failure of 
FRLs. 
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